home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / TPL60N14.ARJ / DSQRT.PAS < prev    next >
Pascal/Delphi Source File  |  1992-05-01  |  8KB  |  266 lines
  1. PROGRAM DSqrt;   { ported from Fortran original 05-01-92 Norbert Juffa }
  2.  
  3. {$A+,B-,D-,E+,F-,G-,I-,L-,N-,O-,R-,S-,V-,X-}
  4.  
  5. USES MachArit;
  6.  
  7. {
  8. C     PROGRAM TO TEST DSQRT
  9. C
  10. C     DATA REQUIRED
  11. C
  12. C        NONE
  13. C
  14. C     SUBPROGRAMS REQUIRED FROM THIS PACKAGE
  15. C
  16. C        MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING
  17. C                 INFORMATION ON THE FLOATING-POINT ARITHMETIC
  18. C                 SYSTEM.  NOTE THAT THE CALL TO MACHAR CAN
  19. C                 BE DELETED PROVIDED THE FOLLOWING SIX
  20. C                 PARAMETERS ARE ASSIGNED THE VALUES INDICATED
  21. C
  22. C                 IBETA  - THE RADIX OF THE FLOATING-POINT SYSTEM
  23. C                 IT     - THE NUMBER OF BASE-IBETA DIGITS IN THE
  24. C                          SIGNIFICAND OF A FLOATING-POINT NUMBER
  25. C                 EPS    - THE SMALLEST POSITIVE FLOATING-POINT
  26. C                          NUMBER SUCH THAT 1.0+EPS .NE. 1.0
  27. C                 EPSNEG - THE SMALLEST POSITIVE FLOATING-POINT
  28. C                          NUMBER SUCH THAT 1.0-EPSNEG .NE. 1.0
  29. C                 XMIN   - THE SMALLEST NON-VANISHING FLOATING-POINT
  30. C                          POWER OF THE RADIX
  31. C                 XMAX   - THE LARGEST FINITE FLOATING-POINT NO.
  32. C
  33. C      RANDL(X) - A FUNCTION SUBPROGRAM RETURNING LOGARITHMICALLY
  34. C                 DISTRIBUTED RANDOM REAL NUMBERS.  IN PARTICULAR,
  35. C                        A * RANDL(DLOG(B/A))
  36. C                 IS LOGARITHMICALLY DISTRIBUTED OVER (A,B)
  37. C
  38. C        REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL
  39. C                 NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1)
  40. C
  41. C
  42. C     STANDARD FORTRAN SUBPROGRAMS REQUIRED
  43. C
  44. C         DABS, DLOG, DMAX1, DFLOAT, DSQRT
  45. C
  46. C
  47. C     LATEST REVISION - AUGUST 2, 1979
  48. C
  49. C     AUTHOR - W. J. CODY
  50. C              ARGONNE NATIONAL LABORATORY
  51. C
  52. C
  53. }
  54.  
  55.  
  56. FUNCTION REN (K: LONGINT): REAL;
  57.  
  58. {
  59.       DOUBLE PRECISION FUNCTION REN(K)
  60. C
  61. C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND
  62. C      HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM,
  63. C      VOL. 8, NO. 10, OCTOBER 1965.
  64. C
  65. C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH
  66. C      FIXED POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS
  67. C      BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST
  68. C      29 BITS.
  69. C
  70. }
  71.  
  72. VAR   J: LONGINT;
  73. CONST IY: LONGINT = 100001;
  74.  
  75. BEGIN
  76.    J  := K;
  77.    IY := IY * 125;
  78.    IY := IY - (IY DIV 2796203) * 2796203;
  79.    REN:= 1.0 * (IY) / 2796203.0e0 * (1.0e0 + 1.0e-6 + 1.0e-12);
  80. END;
  81.  
  82.  
  83. FUNCTION MAX1 (A, B:REAL): REAL;
  84. BEGIN
  85.    IF A > B THEN
  86.       MAX1 := A
  87.    ELSE
  88.       MAX1 := B;
  89. END;
  90.  
  91.  
  92.  
  93. FUNCTION RANDL(X: REAL): REAL;
  94. {
  95. C
  96. C     RETURNS PSEUDO RANDOM NUMBERS LOGARITHMICALLY DISTRIBUTED
  97. C     OVER (1,EXP(X)).  THUS A*RANDL(LN(B/A)) IS LOGARITHMICALLY
  98. C     DISTRIBUTED IN (A,B).
  99. C
  100. C     OTHER SUBROUTINES REQUIRED
  101. C
  102. C        EXP(X) - THE EXPONENTIAL ROUTINE
  103. C
  104. C        REN(K) - A FUNCTION PROGRAM RETURNING RANDOM REAL
  105. C                 NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1).
  106. C                 THE ARGUMENT K IS A DUMMY.
  107. C
  108. C
  109. }
  110.  
  111. CONST K:LONGINT=1;
  112. BEGIN
  113.    RANDL := EXP (X*REN(K));
  114. END;
  115.  
  116.  
  117.  
  118. VAR   I,IBETA,IEXP,IOUT,IRND,IT,J,K1,K2,
  119.       K3,MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD: LONGINT;
  120.       A,AIT,ALBETA,B,BETA,C,EPS,EPSNEG,ONE,
  121.       R6,R7,SQBETA,W,X,XMAX,XMIN,XN,X1,Y,Z,ZERO: REAL;
  122.  
  123. LABEL 100, 110, 120, 150, 160, 210, 220, 230, 240, 300;
  124.  
  125. BEGIN
  126.  
  127.    N := 10000;   { number of trials }
  128.  
  129.    MACHAR (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
  130.            EPS,EPSNEG,XMIN,XMAX);
  131.    PRINTPARAM (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
  132.                EPS,EPSNEG,XMIN,XMAX);
  133.    BETA  := IBETA;
  134.    SQBETA:= SQRT (BETA);
  135.    ALBETA:= LN (BETA);
  136.    AIT   := (IT);
  137.    ONE   := 1;
  138.    ZERO  := 0;
  139.    A     := ONE / SQBETA;
  140.    B     := ONE;
  141.    XN    := N;
  142.  
  143. {-----------------------------------------------------------------}
  144. {     RANDOM ARGUMENT ACCURACY TESTS                              }
  145. {-----------------------------------------------------------------}
  146.  
  147.    FOR J := 1 TO 2 DO BEGIN
  148.       C  := LN (B/A);
  149.       K1 := 0;
  150.       K3 := 0;
  151.       X1 := ZERO;
  152.       R6 := ZERO;
  153.       R7 := ZERO;
  154.  
  155.       FOR I := 1 TO N DO BEGIN
  156.          X := A * RANDL(C);
  157.          Y := X * X;
  158.          Z := SQRT (Y);
  159.          IF X <> ZERO THEN
  160.             W := (Z - X) / X
  161.          ELSE IF Z <> ZERO THEN
  162.             W := ONE;
  163.          IF W > ZERO THEN
  164.             K1 := K1 + 1;
  165.          IF W < ZERO THEN
  166.             K3 := K3 + 1;
  167.          W := ABS (W);
  168.          IF W <= R6 THEN
  169.             GOTO 120;
  170.          R6 := W;
  171.          X1 := X;
  172. 120:     R7 := R7 + W * W;
  173.       END;
  174.  
  175.       K2 := N - K1 - K3;
  176.       R7 := SQRT (R7/XN);
  177.  
  178.       WRITELN;
  179.       WRITELN;
  180.       WRITELN ('TEST OF SQRT(X*X) - X ');
  181.       WRITELN;
  182.       WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
  183.       WRITELN ('(', A, ',', B, ')');
  184.       WRITELN;
  185.       WRITELN ('SQRT (X) WAS LARGER', K1:6, ' TIMES');
  186.       WRITELN ('             AGREED', K2:6, ' TIMES');
  187.       WRITELN ('    AND WAS SMALLER', K3:6, ' TIMES');
  188.       WRITELN;
  189.       WRITELN ('THERE ARE ', IT, ' BASE ', IBETA,
  190.                ' SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER');
  191.       WRITELN;
  192.       W := -999.0;
  193.       IF R6 <> ZERO THEN
  194.          W := LN (ABS(R6))/ALBETA;
  195.       WRITELN ('THE MAXIMUM RELATIVE ERROR OF          ', R6:12,
  196.                ' = ', IBETA, ' **', W:7:2);
  197.       WRITELN ('OCCURED FOR X = ', X1);
  198.       W := MAX1 (AIT+W,ZERO);
  199.       WRITELN;
  200.       WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  201.                ' SIGNIFICANT DIGITS IS        ', W:7:2);
  202.       W := -999;
  203.       IF R7 <> ZERO THEN
  204.          W := LN (ABS(R7))/ALBETA;
  205.       WRITELN;
  206.       WRITELN ('THE ROOT MEAN SQUARE RELATIVE ERROR WAS', R7:12,
  207.                ' = ', IBETA, ' **' , W:7:2);
  208.       W := MAX1 (AIT+W,ZERO);
  209.       WRITELN;
  210.       WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  211.                ' SIGNIFICANT DIGITS IS        ', W:7:2);
  212.       A := ONE;
  213.       B := SQBETA;
  214.    END;
  215.  
  216. {-----------------------------------------------------------------}
  217. {     SPECIAL TESTS                                               }
  218. {-----------------------------------------------------------------}
  219.  
  220.    WRITELN;
  221.    WRITELN;
  222.    WRITELN ('TEST OF SPECIAL ARGUMENTS');
  223.    WRITELN;
  224.    X := XMIN;
  225.    Y := SQRT (X);
  226.    WRITELN ('SQRT (XMIN)    = SQRT (  ', X:18, ') = ', Y:18);
  227.    WRITELN;
  228.    X := ONE - EPSNEG;
  229.    Y := SQRT(X);
  230.    WRITELN ('SQRT(1-EPSNEG) = SQRT (1-', EPSNEG:18, ') = ', Y:18);
  231.    WRITELN;
  232.    X := ONE;
  233.    Y := SQRT(X);
  234.    WRITELN ('SQRT (1.0)     = SQRT (  ', X:18, ') = ', Y:18);
  235.    WRITELN;
  236.    X := ONE + EPS;
  237.    Y := SQRT(X);
  238.    WRITELN ('SQRT (1+EPS)   = SQRT (1+', EPS:18, ') = ', Y:18);
  239.    WRITELN;
  240.    X := XMAX;
  241.    Y := SQRT(X);
  242.    WRITELN ('SQRT (XMAX)    = SQRT (  ', X:18, ') = ', Y:18);
  243.    WRITELN;
  244.  
  245. {-----------------------------------------------------------------}
  246. {     TEST OF ERROR RETURNS                                       }
  247. {-----------------------------------------------------------------}
  248.  
  249.    WRITELN;
  250.    WRITELN;
  251.    WRITELN ('TEST OF ERROR RETURNS');
  252.    WRITELN;
  253.    X := ZERO;
  254.    WRITELN ('SQRT WILL BE CALLED WITH THE ARGUMENT ',  X:15);
  255.    WRITELN ('THIS SHOULD NOT TRIGGER AN ERROR MESSAGE');
  256.    Y := SQRT(X);
  257.    WRITELN ('SQRT RETURNED THE VALUE ', Y:15);
  258.    X := -ONE;
  259.    WRITELN ('SQRT WILL BE CALLED WITH THE ARGUMENT ',  X:15);
  260.    WRITELN ('THIS SHOULD TRIGGER AN ERROR MESSAGE');
  261.    Y := SQRT(X);
  262.    WRITELN ('SQRT RETURNED THE VALUE ', Y:15);
  263.    WRITELN;
  264.    WRITELN ('THIS CONCLUDES THE TESTS');
  265. END. { DSqrt }
  266.